///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file energy_accumulator.cpp
///
/// \brief   implementation for the energy accumulator class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<src/pcps.h>
#include<src/input.h>
#include<src/configuration.h>
#include<src/accumulator.h>
#include<src/energy_accumulator.h>
#include<src/hamiltonian.h>
#include<src/timing.h>

#include<src/restricted_pairing_function.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Constructs the energy accumulator.
///
/// \param[in]      userinp       the user's input options
/// \param[in]      ham           the Hamiltonian
/// \param[in]      co            a pointer to the coefficient object to use
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> pcps::EnergyAccumulator<S>::EnergyAccumulator(const pcps::Input & userinp,
                                                                 const pcps::Hamiltonian<S> & ham,
                                                                 boost::shared_ptr<pcps::CoeffObj<S> > & co)
  : Accumulator<S>(),
    _nas(userinp.nsites()/2),
    _nap(userinp.nparticles()/2),
    _nau(userinp.nsites()/2 - userinp.nparticles()/2),
    _max_ac_length(userinp.max_ac_length()),
    _nch(ham.n_cholesky()),
    _co(co),
    _ham(&ham),
    _wfn_sampling(true),
    _pqpq_guiding(userinp.pqpq_guiding()),
    _print(false),
    _energy(0.0),
    _pqqp_multiplier(0.0),
    _oei_contribution(0.0),
    _pqqp_energy(0.0),
    _pppp_energy(0.0),
    _a_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _b_ratios( (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _a_contraction(ham.n_cholesky(), pcps::zero<S>()),
    _b_contraction(ham.n_cholesky(), pcps::zero<S>()),
    _a_diag_cont(ham.n_cholesky(), pcps::zero<S>()),
    _b_diag_cont(ham.n_cholesky(), pcps::zero<S>()),
    _oei_slice((userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _chol_slice(ham.n_cholesky() * (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _a_oei_slice((userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _a_chol_slice(ham.n_cholesky() * (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _b_oei_slice((userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _b_chol_slice(ham.n_cholesky() * (userinp.nparticles()/2) * ( userinp.nsites()/2 - userinp.nparticles()/2 + 1 ), pcps::zero<S>()),
    _pqpq_tei( userinp.nsites() * userinp.nsites() / 4 ),
    _history(),
    _config(userinp)
{

  // check that we have zero net spin
  if (userinp.z_spin() != 0.0)
    throw pcps::Exception("energy accumulator constructor assumes equal number of alpha and beta electrons");

  // get the pqpq two electron integrals
  ham.get_pqpq_ints(_pqpq_tei);

}

template pcps::EnergyAccumulator<double>::EnergyAccumulator(
 const pcps::Input & userinp,
 const pcps::Hamiltonian<double> & ham,
 boost::shared_ptr<pcps::CoeffObj<double> > & co);
template pcps::EnergyAccumulator<std::complex<double> >::EnergyAccumulator(
 const pcps::Input & userinp,
 const pcps::Hamiltonian<std::complex<double> > & ham,
 boost::shared_ptr<pcps::CoeffObj<std::complex<double> > > & co);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  get the slice of one electron integrals and cholesky vectors for the stored configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::init_slices() {

  pcps::start_timer("init_slice");

  for (int b = 0; b < 2; b++) {

    const std::vector<int> & occ_to_site   = ( b == 0 ? _config.a_occ_to_site()   : _config.b_occ_to_site()   );
    const std::vector<int> & unocc_to_site = ( b == 0 ? _config.a_unocc_to_site() : _config.b_unocc_to_site() );
    std::vector<S> & oei_slice = ( b == 0 ? _a_oei_slice : _b_oei_slice );
    std::vector<S> & chol_slice = ( b == 0 ? _a_chol_slice : _b_chol_slice );

    for (int i = 0; i < occ_to_site.size(); i++) {

      const int p = occ_to_site.at(i);

      for (int j = -1; j < int(unocc_to_site.size()); j++) {

        const int q = ( j < 0 ? p : unocc_to_site.at(j) );

        // get the compound index
        const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*q+b);

        // get the modified one electron integral
        oei_slice.at(cpq) = _ham->get_mei(p,q);
        //oei_slice.at(cpq) = _ham->get_oei(p,q);

        // get the necessary elements of the cholesky vectors
        pcps::xcopy(_nch, &_ham->cholesky_vecs_t().at(_nch*(p*_nas+q)), 1, &chol_slice.at(_nch*cpq), 1);

      }

    }

  }

  pcps::stop_timer("init_slice");

}

template void pcps::EnergyAccumulator<double>::init_slices();
template void pcps::EnergyAccumulator<std::complex<double> >::init_slices();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  update the internal configuration and associated slices
///
/// \param[in]      config        the new configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::update_slices(const pcps::Configuration<S> & config) {

  // get the difference between the new and stored configurations
  pcps::ConfigDiff<S> diff(config, _config);

  // if the configurations are the same, then we don't need to update anything
  if ( diff.newly_empty().size() == 0 && diff.newly_full().size() == 0 ) return;

  // make sure the difference is what we expected
  assert( diff.newly_empty().size() == 1 );
  assert( diff.newly_full().size() == 1 );
  assert( diff.newly_full().at(0) % 2 == diff.newly_empty().at(0) % 2 );

  // update the internal configuration
  _config.update(diff);

  //this->init_slices();
  //return;

  // get the spin of the electron move
  const int b = diff.newly_full().at(0) % 2;

  // get the occupancy/unnocupancy maps and the slice vectors of the approproate spin
  const std::vector<int> & occ_to_site   = ( b == 0 ? _config.a_occ_to_site()   : _config.b_occ_to_site()   );
  const std::vector<int> & unocc_to_site = ( b == 0 ? _config.a_unocc_to_site() : _config.b_unocc_to_site() );
  std::vector<S> & oei_slice = ( b == 0 ? _a_oei_slice : _b_oei_slice );
  std::vector<S> & chol_slice = ( b == 0 ? _a_chol_slice : _b_chol_slice );

  // update the affected slices for each occupied orbital
  for (int i = 0; i < occ_to_site.size(); i++) {

    // get the excitation's indices
    const int p = occ_to_site.at(i);
    const int q = diff.newly_empty().at(0) / 2;

    // get the compound index
    const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*q+b);

    // get the modified one electron integral
    oei_slice.at(cpq) = _ham->get_mei(p,q);
    //oei_slice.at(cpq) = _ham->get_oei(p,q);

    // get the necessary elements of the cholesky vectors
    pcps::xcopy(_nch, &_ham->cholesky_vecs_t().at(_nch*(p*_nas+q)), 1, &chol_slice.at(_nch*cpq), 1);

  }

  // update the affected slices for each unoccupied orbital and the diagonal term
  for (int j = -1; j < int(unocc_to_site.size()); j++) {

    // get the excitation's indices
    const int p = diff.newly_full().at(0) / 2;
    const int q = ( j < 0 ? p : unocc_to_site.at(j) );

    // get the compound index
    const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*q+b);

    // get the modified one electron integral
    oei_slice.at(cpq) = _ham->get_mei(p,q);
    //oei_slice.at(cpq) = _ham->get_oei(p,q);

    // get the necessary elements of the cholesky vectors
    pcps::xcopy(_nch, &_ham->cholesky_vecs_t().at(_nch*(p*_nas+q)), 1, &chol_slice.at(_nch*cpq), 1);

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  get the slice of one electron integrals and cholesky vectors for the given configuration
///
/// \param[in]      userinp       the user's input options
/// \param[in]      config        the configuration that determines the slice
/// \param[in]      b             slice is for alpha or beta,  0-alpha,  1-beta
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::get_ham_slice(const pcps::Input & userinp,
                                                                  const pcps::Configuration<S> & config,
                                                                  const int b) {

  pcps::start_timer("get_ham_slice");

  // check that we have zero net spin
  if (userinp.z_spin() != 0.0)
    throw pcps::Exception("energy accumulator cholesky slice assumes equal number of alpha and beta electrons");

  // get a reference to the fock configuration vector
  const std::vector<int> & fc = config.fock_config();

  // initialize the row counter
  int p_count = 0;

  // index for checking that we are looping correctly
  int last_pq = -1;

  // loop over the p operator
  for (int p = 0; p < _nas; p++) {

    // skip excitations that destroy the bra
    if ( fc[2*p+b] == 0 ) continue;

    // initialize the column counter
    int q_count = 0;

    // loop over the q operator
    for (int q = 0; q < _nas; q++) {

      // skip excitations that destroy the bra
      if ( q != p && fc[2*q+b] == 1 ) continue;

      // get the compound index
      const int cpq = (_nau+1) * p_count + q_count;

      // get the modified one electron integral
      _oei_slice.at(cpq) = _ham->get_mei(p,q);
      //_oei_slice.at(cpq) = _ham->get_oei(p,q);

      // get the necessary elements of the cholesky vectors
      pcps::xcopy(_nch, &_ham->cholesky_vecs_t().at(_nch*(p*_nas+q)), 1, &_chol_slice.at(_nch*cpq), 1);

      // increment the ratio storage column counter
      q_count++;

      // check that we are looping correctly
      assert( cpq == last_pq + 1 );
      last_pq = cpq;

    }

    // increment the ratio storage row counter
    p_count++;

  }

  // check that the correct number of elements were sliced
  assert( last_pq + 1 == (_nau+1)*_nap );

  pcps::stop_timer("get_ham_slice");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Get the alpha or beta contribution to the local energy.
///
/// \param[in]      userinp       the user's input options
/// \param[in]      w             the weight to apply to this configuration's contribution,
///                               including the square of the wavefunction to guiding ratio
/// \param[in]      config        the configuration to accumulate data for
/// \param[in]      b             0 -- alpha     1 -- beta
/// \param[out]     contraction   on exit, the contraction of the alpha/beta ratio matrix
///                               with the cholesky vectors
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::calculate_a(const pcps::Input & userinp,
                                                                const double w,
                                                                const pcps::Configuration<S> & config,
                                                                const int b,
                                                                const std::vector<S> & ratios,
                                                                std::vector<S> & contraction,
                                                                std::vector<S> & diag_cont,
                                                                const std::vector<S> & opp_spin_ratios)
{

  // check that we have zero net spin
  if (userinp.z_spin() != 0.0)
    throw pcps::Exception("energy accumulator currently assumes equal number of alpha and beta electrons");

  //// get the slices of the one electron integrals and cholesky vectors that we need
  //this->get_ham_slice(userinp, config, b);

  //const std::vector<S> & oei_slice = _oei_slice;
  //const std::vector<S> & chol_slice = _chol_slice;
  const std::vector<S> & oei_slice = ( b == 0 ? _a_oei_slice : _b_oei_slice );
  const std::vector<S> & chol_slice = ( b == 0 ? _a_chol_slice : _b_chol_slice );

  if (_print) {

    const std::vector<int> & occ_to_site   = ( b == 0 ? _config.a_occ_to_site()   : _config.b_occ_to_site()   );
    const std::vector<int> & unocc_to_site = ( b == 0 ? _config.a_unocc_to_site() : _config.b_unocc_to_site() );

    const std::vector<int> & os_site_to_occ   = ( b == 0 ? _config.b_site_to_occ()   : _config.a_site_to_occ()   );
    const std::vector<int> & os_site_to_unocc = ( b == 0 ? _config.b_site_to_unocc() : _config.a_site_to_unocc() );

    // if printing, get the contribution to the pqqp energy, p != q
    if (!_wfn_sampling)
      for (int i = 0; i < occ_to_site.size(); i++) {
        const int p = occ_to_site[i];
        for (int j = 0; j < unocc_to_site.size(); j++) {
          const int q = unocc_to_site[j];
          const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*q+b);
          // aa contribution
          for (int u = 0; u < _nch; u++)
            _pqqp_energy += 0.5 * pcps::square_norm( _ham->cholesky_vecs_t().at(_nch*(p*_nas+q) + u) * ratios.at(cpq) );
            //_pqqp_energy += pcps::square_norm( chol_slice.at( (_nau+1)*_nap * u + cpq ) * ratios.at(cpq) );
          // ab contribution
          const int a = (b+1)%2;
          int os_cpq = -1;
          if ( os_site_to_occ[p] >= 0 && os_site_to_unocc[q] >= 0 )
            os_cpq = _config.occ_unocc_cmpd(2*p+a, 2*q+a);
          if ( os_site_to_occ[q] >= 0 && os_site_to_unocc[p] >= 0 )
            os_cpq = _config.occ_unocc_cmpd(2*q+a, 2*p+a);
          if (os_cpq >= 0)
            for (int u = 0; u < _nch; u++)
              _pqqp_energy += 0.5 * pcps::square_norm( _ham->cholesky_vecs_t().at(_nch*(p*_nas+q) + u) )
                                  * pcps::real( ratios.at(cpq) * opp_spin_ratios.at(os_cpq) );
        }
      }

    // if printing, get the contribution to the left pqqp energy, p != q, <Psi|n><n|pqqp|Psi>
    if (_wfn_sampling)
      for (int i = 0; i < occ_to_site.size(); i++) {
        const int p = occ_to_site[i];
        for (int j = 0; j < unocc_to_site.size(); j++) {
          const int q = unocc_to_site[j];
          const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*q+b);
          for (int u = 0; u < _nch; u++)
            _pqqp_energy += 0.5 * pcps::square_norm( _ham->cholesky_vecs_t().at(_nch*(p*_nas+q) + u) );
        }
      }

    // if printing, get the contribution to the pppp energy
    if (!_wfn_sampling)
      for (int i = 0; i < occ_to_site.size(); i++) {
        const int p = occ_to_site[i];
        const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*p+b);
        for (int u = 0; u < _nch; u++)
          _pppp_energy += pcps::square_norm( _ham->cholesky_vecs_t().at(_nch*(p*_nas+p) + u) * ratios.at(cpq) );
          //_pppp_energy += pcps::square_norm( chol_slice.at( (_nau+1)*_nap * u + cpq ) * ratios.at(cpq) );
      }

  }

  pcps::start_timer("big contraction");

  // get the contraction of the ratios with the cholesky vectors
  pcps::xgemm('N', 'N', _nch, 1, (_nau+1)*_nap, pcps::unity<S>(),
              &chol_slice.at(0), _nch, &ratios.at(0), (_nau+1)*_nap,
              pcps::zero<S>(), &contraction.at(0), _nch);

  pcps::stop_timer("big contraction");

  pcps::start_timer("diag contraction");

  // get the diagonal contraction
  if (true) {

    std::fill(diag_cont.begin(), diag_cont.end(), pcps::zero<S>());

    const std::vector<int> & occ_to_site   = ( b == 0 ? _config.a_occ_to_site()   : _config.b_occ_to_site()   );

    for (int i = 0; i < occ_to_site.size(); i++) {
      const int p = occ_to_site[i];
      const int cpq = _config.occ_unocc_cmpd(2*p+b, 2*p+b);
      pcps::xaxpy(_nch, ratios.at(cpq), &chol_slice.at(_nch*cpq), 1, &diag_cont.at(0), 1);
    }

  // get the diagonal contraction
  } else {

    std::fill(diag_cont.begin(), diag_cont.end(), pcps::zero<S>());

    // get a reference to the fock configuration vector
    const std::vector<int> & fc = config.fock_config();

    // initialize the row counter
    int p_count = 0;

    // index for checking that we are looping correctly
    int last_pq = -1;

    // loop over the p operator
    for (int p = 0; p < _nas; p++) {

      // skip excitations that destroy the bra
      if ( fc[2*p+b] == 0 ) continue;

      // initialize the column counter
      int q_count = 0;

      // loop over the q operator
      for (int q = 0; q < _nas; q++) {

        // skip excitations that destroy the bra
        if ( q != p && fc[2*q+b] == 1 ) continue;

        // get the compound index
        const int cpq = (_nau+1) * p_count + q_count;

        // add this site's diagonal contribution
        if ( p == q )
          pcps::xaxpy(_nch, ratios.at(cpq), &chol_slice.at(_nch*cpq), 1, &diag_cont.at(0), 1);

        // increment the ratio storage column counter
        q_count++;

        // check that we are looping correctly
        assert( cpq == last_pq + 1 );
        last_pq = cpq;

      }

      // increment the ratio storage row counter
      p_count++;

    }

    // check that the correct number of elements were looped over
    assert( last_pq + 1 == (_nau+1)*_nap );

  }

  pcps::stop_timer("diag contraction");

  pcps::start_timer("other calculate_a");

  // wavefunction sampling
  if (_wfn_sampling) {

    // get the contribution from the one electron integrals
    const double x = w * pcps::real( pcps::xdot((_nau+1)*_nap, &oei_slice.at(0), 1, &ratios.at(0), 1) );
    _oei_contribution += x;
    _energy += x;

    // get the contribution to the pqqp multiplier
    const std::vector<int> & occ_to_site   = ( b == 0 ? _config.a_occ_to_site()   : _config.b_occ_to_site()   );
    const std::vector<int> & unocc_to_site = ( b == 0 ? _config.a_unocc_to_site() : _config.b_unocc_to_site() );
    for (int i = 0; i < occ_to_site.size(); i++) {
      const int p = occ_to_site[i];
      for (int j = 0; j < unocc_to_site.size(); j++) {
        const int q = unocc_to_site[j];
        _pqqp_multiplier += 0.5 * w * pcps::real(_pqpq_tei[_nas*p+q]);
      }
    }

  // rdm sampling
  } else {

    // remove the diagonal contribution to the contraction
    pcps::xaxpy(_nch, -pcps::unity<S>(), &diag_cont.at(0), 1, &contraction.at(0), 1);

  }

  pcps::stop_timer("other calculate_a");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Accumulate the local energy (note this is the sandwiched cholesky decomposition version
///        and not the typical VMC local energy).
///
/// \param[in]      userinp  the user's input options
/// \param[in]      w        the weight to apply to this configuration's contribution
/// \param[in]      config   the configuration to accumulate data for
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::child_accumulate(const pcps::Input & userinp,
                                                                     const double w,
                                                                     const pcps::Configuration<S> & config)
{

  // zero some of the to-be-printed values
  _oei_contribution = 0.0;
  _pqqp_energy = 0.0;
  _pppp_energy = 0.0;

  // get the old energy
  const double old_energy = _energy;

  // get the single excitation ratios
  _co->single_excitation_ratios(config, _a_ratios, _b_ratios);

  pcps::start_timer("update_slices");

  // update the slices
  if (_print) {
    _config = config;
    this->init_slices();
  } else
    this->update_slices(config);

  pcps::stop_timer("update_slices");

  //pcps::start_timer("calculate_a");

  // for alpha excitations, compute cholesky/ratios contractions and add the one electron energy contribution
  this->calculate_a(userinp, w, config, 0, _a_ratios, _a_contraction, _a_diag_cont, _b_ratios);

  // for beta excitations, compute cholesky/ratios contractions and add the one electron energy contribution
  this->calculate_a(userinp, w, config, 1, _b_ratios, _b_contraction, _b_diag_cont, _a_ratios);

  //pcps::stop_timer("calculate_a");

  pcps::start_timer("contract contractions");

  // wavefunction sampling
  if (_wfn_sampling) {

    // add the two electron energy contribution for the half-diagonal two electron integrals
    const double half_diag_cont = 0.5 * w * pcps::real(   2.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0), 1, &_a_contraction.at(0), 1)
                                                        - 1.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0), 1, &_a_diag_cont.at(0),   1)
                                                        + 2.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0), 1, &_b_contraction.at(0), 1)
                                                        - 1.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0), 1, &_b_diag_cont.at(0),   1)
                                                        + 2.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0), 1, &_b_contraction.at(0), 1)
                                                        + 2.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0), 1, &_a_contraction.at(0), 1)
                                                        - 2.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0), 1, &_b_diag_cont.at(0),   1) );

    _energy += half_diag_cont;

    if (_print) {
      std::cout << boost::format("  %24.12f") % _co->get_value().get();
      std::cout << boost::format("  %24.12f") % ( _oei_contribution * _co->get_value().get() );
      std::cout << boost::format("  %24.12f") % ( half_diag_cont * _co->get_value().get() * _co->get_value().get() );
      std::cout << boost::format("  %24.12f") % ( _pqqp_energy * _co->get_value().get() * _co->get_value().get() );
    }

  // rdm sampling
  } else {

    // add the two electron energy contribution for the off-diagonal two electron integrals
    const double x = 0.5 * pcps::real( (         pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_a_contraction.at(0), 1)
                                         +       pcps::xdotc(_nch, &_b_contraction.at(0), 1, &_b_contraction.at(0), 1)
                                         + 2.0 * pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_b_contraction.at(0), 1) ) / _co->get_value().get() );

//    const pcps::RestPairingFunc<S> * const rpf = dynamic_cast<pcps::RestPairingFunc<S> *>( _co.get() );
//    const double cor_val = pcps::real( rpf->correct_value(config, 'e') );
//    //std::cout << boost::format("   local energy  %20.12f    correct  %20.12f")
//    std::cout << boost::format("   energy error  %20.12f")
//                 % ( x - cor_val )
//              << std::endl;
//    //assert( std::abs( x - cor_val ) < 1.0e-9 );
//
//    _energy += cor_val * w;
//
//    std::cout << boost::format("local energy    %20.12f") % (x/w) << std::endl;

    _energy += x * w * _pqqp_multiplier;

    if (_print) {
      const S init_val = ( dynamic_cast<pcps::RestPairingFunc<S> *>(&*_co)  )->init_value().get();
      std::cout << boost::format("  %24.12f") % ( _co->get_value().get() * init_val * init_val );
      std::cout << boost::format("  %24.12f") % (                      x * init_val * init_val * _co->get_value().get() ); // off-diagonal tei energy
      std::cout << boost::format("  %24.12f") % (           _pqqp_energy * init_val * init_val );
      std::cout << boost::format("  %24.12f") % (           _pppp_energy * init_val * init_val );
      const double xaa = 0.5 * pcps::real( pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_a_contraction.at(0), 1) * init_val * init_val );
      const double xbb = 0.5 * pcps::real( pcps::xdotc(_nch, &_b_contraction.at(0), 1, &_b_contraction.at(0), 1) * init_val * init_val );
      const double xab = 1.0 * pcps::real( pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_b_contraction.at(0), 1) * init_val * init_val );
      std::cout << boost::format("  %24.12f") % ( xaa );
      std::cout << boost::format("  %24.12f") % ( xbb );
      std::cout << boost::format("  %24.12f") % ( xab );
    }

//    _energy += 0.5 * w * pcps::real( (   1.0 * pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_a_contraction.at(0), 1)
//                                       + 1.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0),   1, &_a_diag_cont.at(0),   1)
//                                       - 2.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0),   1, &_a_contraction.at(0), 1)
//
//                                       + 1.0 * pcps::xdotc(_nch, &_b_contraction.at(0), 1, &_b_contraction.at(0), 1)
//                                       + 1.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0),   1, &_b_diag_cont.at(0),   1)
//                                       - 2.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0),   1, &_b_contraction.at(0), 1)
//
//                                       + 1.0 * pcps::xdotc(_nch, &_a_contraction.at(0), 1, &_b_contraction.at(0), 1)
//                                       + 1.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0),   1, &_b_diag_cont.at(0),   1)
//                                       - 2.0 * pcps::xdotc(_nch, &_a_diag_cont.at(0),   1, &_b_contraction.at(0), 1)
//
//                                       + 1.0 * pcps::xdotc(_nch, &_b_contraction.at(0), 1, &_a_contraction.at(0), 1)
//                                       + 1.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0),   1, &_a_diag_cont.at(0),   1)
//                                       - 2.0 * pcps::xdotc(_nch, &_b_diag_cont.at(0),   1, &_a_contraction.at(0), 1) ) / _co->get_value().get() );

  }

  pcps::stop_timer("contract contractions");

  // update the history of samples
  for (int i = pcps::round(w); i > 0; i--)
    _history.push_back( ( _energy - old_energy ) / w );

}

template void pcps::EnergyAccumulator<double>::child_accumulate(const pcps::Input &,
                                                                const double,
                                                                const pcps::Configuration<double> &);
template void pcps::EnergyAccumulator<std::complex<double> >::child_accumulate(const pcps::Input &,
                                                                               const double,
                                                                               const pcps::Configuration<std::complex<double> > &);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Prints a header for the detailed print information
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::print_header() {

  if (_wfn_sampling) {
    std::cout << boost::format("  %24s") % "<n|Psi>";
    std::cout << boost::format("  %24s") % "<n|H_1|Psi>";
    std::cout << boost::format("  %24s") % "pprs \"half-diagonal\"";
    std::cout << boost::format("  %24s") % "<Psi|n><n|pqqp|Psi>";
  } else {
      std::cout << boost::format("  %24.12f") % "2rdm weight";
      std::cout << boost::format("  %24.12f") % "pqrs \"off-diagonal\"";
      std::cout << boost::format("  %24.12f") % "pqqp energy, p!=q";
      std::cout << boost::format("  %24.12f") % "pppp energy";
      std::cout << boost::format("  %24.12f") % "pqrs aa \"off-diagonal\"";
      std::cout << boost::format("  %24.12f") % "pqrs bb \"off-diagonal\"";
      std::cout << boost::format("  %24.12f") % "pqrs ab \"off-diagonal\"";
  }

}

template void pcps::EnergyAccumulator<double>::print_header();
template void pcps::EnergyAccumulator<std::complex<double> >::print_header();

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Prints some statistics about the most recent sample.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::EnergyAccumulator<S>::print_statistics() {

  
  // compute the mean and variance
  const double mean = std::accumulate(_history.begin(), _history.end(), 0.0) / _history.size();

  // compute the variance
  const double var = std::accumulate(_history.begin(), _history.end(), 0.0, VarBinOp(mean)) / _history.size();

  std::cout << boost::format("energy accumulation statistics on process %i:") % pcps::mpi::rank() << std::endl;
  std::cout << std::endl;
  std::cout << boost::format("   acceptance ratio = %20.12f") % (this->_naccept / this->_total_weight) << std::endl;
  std::cout << boost::format("               mean = %20.12f") % mean << std::endl;
  std::cout << boost::format("           variance = %20.12f") % var << std::endl;
  std::cout << boost::format("     sqrt(variance) = %20.12f") % std::sqrt(var) << std::endl;
  std::cout << std::endl;

  // compute the autocorrelation for short times
  std::vector<double> ac(_max_ac_length, 0.0);
  for (int i = 1; i <= ac.size(); i++) {
    ac.at(i-1) = std::accumulate(_history.begin()+i, _history.end(), 0.0, AutoCorrBinOp(mean, i)) / ( _history.size() - i ) / var;
    std::cout << boost::format("autocorrelation %3i = %20.12f") % i % ac.at(i-1) << std::endl;
  }
  std::cout << std::endl;

}

template void pcps::EnergyAccumulator<double>::print_statistics();
template void pcps::EnergyAccumulator<std::complex<double> >::print_statistics();
